#!/usr/bin/env julia

using Printf
using Statistics
using LinearAlgebra
using Random

Random.seed!(42)

# Debug variance scaling issue
function debug_variance_scaling()
    println("DEBUGGING VARIANCE SCALING ISSUE")
    println("="^50)

    # Create ensemble and static components
    state_size = 50  # Smaller for debugging
    ensemble_size = 40

    # Generate synthetic ensemble perturbations with higher variance
    ensemble_perts = 2.0 * randn(state_size, ensemble_size)  # Scale up ensemble variance
    ensemble_cov = (ensemble_perts * ensemble_perts') / (ensemble_size - 1)

    # Create static covariance (diagonal for simplicity)
    static_variances = 1.0 .+ 0.5 * rand(state_size)
    static_cov = Diagonal(static_variances)

    # Compute variances
    ensemble_variance = sum(diag(ensemble_cov))
    static_variance = sum(diag(static_cov))

    println("Ensemble total variance: $(ensemble_variance)")
    println("Static total variance: $(static_variance)")
    println("Ensemble/Static ratio: $(ensemble_variance / static_variance)")
    println()

    # Test hybrid combinations
    hybrid_coeffs = [0.0, 0.25, 0.5, 0.75, 1.0]
    variance_results = Dict{Float64, Any}()

    println("Hybrid coefficient -> Total variance")
    println("-"^40)

    for coeff in hybrid_coeffs
        # Construct hybrid covariance
        hybrid_cov = coeff * ensemble_cov + (1 - coeff) * static_cov

        # Compute eigenvalue decomposition
        eigenvals = eigvals(Matrix(hybrid_cov))

        # Check variance properties
        total_variance = sum(eigenvals)
        variance_results[coeff] = total_variance

        println("$coeff -> $total_variance")
    end

    println()
    println("Checking monotonicity:")
    prev_var = nothing
    monotonic = true
    for coeff in sort(collect(keys(variance_results)))
        curr_var = variance_results[coeff]
        if prev_var !== nothing && curr_var < prev_var
            println("  Non-monotonic at $coeff: $curr_var < $prev_var")
            monotonic = false
        end
        prev_var = curr_var
    end

    if monotonic
        println("  ✓ Scaling is monotonic")
    else
        println("  ✗ Scaling is NOT monotonic")
    end

    println()
    println("DIAGNOSTIC: Ensemble covariance properties:")
    println("  Trace: $(tr(ensemble_cov))")
    println("  Frobenius norm: $(norm(ensemble_cov))")
    println("  Condition number: $(cond(ensemble_cov))")
    println("  Minimum eigenvalue: $(minimum(eigvals(Matrix(ensemble_cov))))")
    println("  Maximum eigenvalue: $(maximum(eigvals(Matrix(ensemble_cov))))")

    println()
    println("DIAGNOSTIC: Static covariance properties:")
    println("  Trace: $(tr(static_cov))")
    println("  Frobenius norm: $(norm(static_cov))")
    println("  Condition number: $(cond(static_cov))")
    println("  Minimum eigenvalue: $(minimum(eigvals(Matrix(static_cov))))")
    println("  Maximum eigenvalue: $(maximum(eigvals(Matrix(static_cov))))")
end

debug_variance_scaling()